Library


library(AICcmodavg)
library(ape)
library(caper)
library(geiger)
library(kableExtra)
library(MuMIn)
library(nlme)
library(phytools)
library(plotrix)
library(shape)
library(car)
library(magick)

Session information


R.version
                 _                           
  platform       x86_64-apple-darwin17.0     
  arch           x86_64                      
  os             darwin17.0                  
  system         x86_64, darwin17.0          
  status                                     
  major          4                           
  minor          2.2                         
  year           2022                        
  month          10                          
  day            31                          
  svn rev        83211                       
  language       R                           
  version.string R version 4.2.2 (2022-10-31)
  nickname       Innocent and Trusting
sessionInfo()
  R version 4.2.2 (2022-10-31)
  Platform: x86_64-apple-darwin17.0 (64-bit)
  Running under: macOS Catalina 10.15.7
  
  Matrix products: default
  BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
  LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
  
  locale:
  [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
  
  attached base packages:
  [1] stats     graphics  grDevices utils     datasets  methods   base     
  
  other attached packages:
   [1] magick_2.7.4     car_3.1-1        carData_3.0-5    shape_1.4.6     
   [5] plotrix_3.8-2    phytools_1.5-12  maps_3.4.1       nlme_3.1-162    
   [9] MuMIn_1.47.5     kableExtra_1.3.4 geiger_2.0.10    caper_1.0.1     
  [13] mvtnorm_1.1-3    MASS_7.3-58.3    ape_5.7-1        AICcmodavg_2.3-2
  [17] rmarkdown_2.20.1
  
  loaded via a namespace (and not attached):
   [1] Rcpp_1.0.10             subplex_1.8             svglite_2.1.1          
   [4] lattice_0.20-45         digest_0.6.31           unmarked_1.2.5         
   [7] foreach_1.5.2           R6_2.5.1                plyr_1.8.8             
  [10] stats4_4.2.2            evaluate_0.20           coda_0.19-4            
  [13] httr_1.4.5              rlang_1.1.0             rstudioapi_0.14        
  [16] jquerylib_0.1.4         phangorn_2.11.1         Matrix_1.5-3           
  [19] combinat_0.0-8          splines_4.2.2           webshot_0.5.4          
  [22] stringr_1.5.0           munsell_0.5.0           igraph_1.4.1           
  [25] compiler_4.2.2          numDeriv_2016.8-1.1     xfun_0.37              
  [28] systemfonts_1.0.4       pkgconfig_2.0.3         mnormt_2.1.1           
  [31] optimParallel_1.0-2     htmltools_0.5.4         expm_0.999-7           
  [34] quadprog_1.5-8          codetools_0.2-19        viridisLite_0.4.1      
  [37] grid_4.2.2              lifecycle_1.0.3         jsonlite_1.8.4         
  [40] xtable_1.8-4            magrittr_2.0.3          scales_1.2.1           
  [43] stringi_1.7.12          cli_3.6.0               cachem_1.0.7           
  [46] pbapply_1.7-0           scatterplot3d_0.3-43    doParallel_1.0.17      
  [49] xml2_1.3.3              bslib_0.4.2             generics_0.1.3         
  [52] fastmatch_1.1-3         deSolve_1.35            iterators_1.0.14       
  [55] tools_4.2.2             glue_1.6.2              abind_1.4-5            
  [58] parallel_4.2.2          fastmap_1.1.1           survival_3.5-5         
  [61] yaml_2.3.7              colorspace_2.1-0        rvest_1.0.3            
  [64] VGAM_1.1-8              clusterGeneration_1.3.7 knitr_1.42             
  [67] sass_0.4.5

Consider only sit-and-wait and widely-foraging lizards

Almost a year ago, my mentors and I submitted a paper to Proceedings of the Royal Society B. In this paper, we investigated the relationship between foraging behavior and reproducitve effort in lizards (see here). One of the most important comments I got from a reviwer was related to the scatter in the data. The reviwer stated that the variance of the data was very large, preventing us from making robust conclusions. In response to this concern, I would like to show here one way to deal with this issue:

As observed in figure 3a-b of the paper, the majority of the species are clustered at small sizes along the x-axis, while a few others are grouped at large sizes and seem to be “extreme cases”. Some of the latter cases are “mixed-foraging” species that can behave either as a “widely-foraging” or “sit-and-wait” species. If I remove the “mixed” category and make the species behave randomly as “widely foragers” or “sit-and-wait forager”, I could recalculate the significance of the interaction to determine whether the results hold or not. If I do so multiple time (>1000), I can calculate the proportion of time that “widely-foraging” species have greater reproductive output than “sit-and-wait” species, as indicated by the original results. Let’s do it:


## Considering only sit-and-wait and widely-foraging lizards


set.seed(94)

options(scipen = 999)

active <- dat_full[dat_full$foraging.mode == "widely foraging",  ]
ambush <- dat_full[dat_full$foraging.mode == "sit and wait",  ]
mixed <- dat_full[dat_full$foraging.mode == "mixed",  ]

savep <- rep(NA, 1000)
savei <- rep(NA, 1000)

for(i in 1:1000){
    
    randomizing <- sample(rep(c("widely foraging", "sit and wait"), each = length(mixed$foraging.mode)), size = length(mixed$foraging.mode))
    mixed$foraging.mode <- randomizing
    merg <- rbind(active, ambush, mixed)

    ## model
    mod <- gls(rep_out ~ M_females*foraging.mode, correlation = corPagel(value = 0, phy = ctree_full, form = ~Binomial), data = merg, method = "ML")

    ## Confidence intervals

    ## sit and wait

    sw <- merg[merg$foraging.mode == "sit and wait", ]
    check_sw <- name.check(tree, sw)
    rm_phy_sw <- check_sw$tree_not_data
    rm_dat_sw <- check_sw$data_not_tree
    ctree_sw <- drop.tip(tree, rm_phy_sw)

    SSX.sw <- sum(round((sw$M_females - mean(sw$M_females))^2), 2)
    s2.sw <- var(sw$rep_out)
    n.sw <- length(sw$rep_out)
    x.sw <- seq(min(sw$M_females), max(sw$M_females), 0.5)
    m.x.sw <- mean(round(sw$M_females, 1))
    se.sw <- sqrt(s2.sw*((1/n.sw) + (((x.sw - m.x.sw)^2)/SSX.sw)))
    is.sw <- qt(0.975, df = n.sw - 2)
    ii.sw <- qt(0.025, df = n.sw - 2)
    ic.s.sw <- se.sw*is.sw
    ic.i.sw <- se.sw*ii.sw
    upper.i.sw <- ((coef(mod)[1] + coef(mod)[3]) + (coef(mod)[2] + coef(mod)[4])*x.sw) + ic.s.sw
    lower.i.sw <- ((coef(mod)[1] + coef(mod)[3]) + (coef(mod)[2] + coef(mod)[4])*x.sw) + ic.i.sw

    ## widely foraging

    wf <- merg[merg$foraging.mode == "widely foraging", ]
    check_wf <- name.check(tree, wf)
    rm_phy_wf <- check_wf$tree_not_data
    rm_dat_wf <- check_wf$data_not_tree
    ctree_wf <- drop.tip(tree, rm_phy_wf)

    SSX.wf <- sum(round((wf$M_females - mean(wf$M_females))^2), 2)
    s2.wf <- var(wf$rep_out)
    n.wf <- length(wf$rep_out)
    x.wf <- seq(min(wf$M_females), max(wf$M_females), 0.5)
    m.x.wf <- mean(round(wf$M_females,1))
    se.wf <- sqrt(s2.wf*((1/n.wf) + (((x.wf - m.x.wf)^2)/SSX.wf)))
    is.wf <- qt(0.975, df = n.wf - 2)
    ii.wf <- qt(0.025, df = n.wf - 2)
    ic.s.wf <- se.wf*is.wf
    ic.i.wf <- se.wf*ii.wf
    upper.i <- (coef(mod)[1] + coef(mod)[2]*x.wf) + ic.s.wf
    lower.i <- (coef(mod)[1] + coef(mod)[2]*x.wf) + ic.i.wf

    png(paste("figure", i, ".png", sep = ""), width = 7, height = 7, units = "in", res = 300)
    plot(merg$rep_out ~ merg$M_females, type = "p", xlab = " ", ylab = "Hatchling mass x clutch size (g)", cex.lab = 1.2, pch = 21, bg = c("purple", "skyblue")[as.numeric(merg$foraging.mode)], col = c("purple", "skyblue")[as.numeric(merg$foraging.mode)], las = 1)
    polygon(c(rev(x.sw), x.sw), c(rev(lower.i.sw), upper.i.sw), border = FALSE, col = rgb(0, 0, 1, alpha = 0.2))
    lines(x = x.sw, y = ((coef(mod)[1] + coef(mod)[3]) + (coef(mod)[2] + coef(mod)[4])*x.sw), col = "skyblue", lwd = 2)

    polygon(c(rev(x.wf), x.wf), c(rev(lower.i), upper.i), border = FALSE, col = rgb(0.8, 0, 0.7, alpha = 0.3))
    lines(x = x.wf, y = coef(mod)[1] + coef(mod)[2]*x.wf, col = "purple", lwd = 2)

    text(x = 7.7, y = 30, substitute(paste(italic("Interaction \n "))), family = "serif")
    text(x = 7.7, y = 30, round(coef(mod)[4], 3), family = "serif")
    text(x = 7.7, y = 27, paste("p = ", round(anova(mod)[4, 3], 3), sep = " "), family = "serif")

    mtext(expression(hat(M) (g)), side = 1, line = 4, cex = 1)
    legend("topleft", legend = c("sit-and-wait foraging", "widely foraging"), lty = 1, bty = "n", pch = c(16, 16), bg = c("black", "black"), col = c("skyblue", "purple"), cex = 1)
    
    dev.off()
    
    
    savep[i] <- round(anova(mod)[4, 3], 3)
    savei[i] <- round(coef(mod)[4], 3)

}
    



## Estimate proportion of time that the p-value of the interaction was significant

mean(savep < 0.05)
  [1] "0.574"

## Mean affect size of the interaction

mean(savei)
  [1] "-0.259"

## list file names and read in

imgs <- list.files(path = "/Users/dpadil10/Dropbox (ASU)/gif/gif/", full.names = TRUE)
img_list <- lapply(imgs, image_read)

## join the images together
img_joined <- image_join(img_list)

## animate at 2 frames per second
img_animated <- image_animate(img_joined, fps = 2)

img_animated